2021: Gulf of Maine’s Warmest Year on Record

Characteristics of a record setting year for ocean temperatures

true
Updated on: 2022-01-21
hide
# Set ggplot theme for figures
theme_set(theme_bw() + 
            theme(
              # Titles
              plot.title = element_text(hjust = 0, face = "bold", size = 14),
              plot.subtitle = element_text(size = 9),
              plot.caption = element_text(size = 7.2, margin = margin(t = 20), color = "gray40"),
              legend.title = element_text(size = 9),
              legend.text = element_text(size = 7.5),
              # Axes
              axis.line.y = element_line(),
              axis.ticks.y = element_line(), 
              axis.line.x = element_line(),
              axis.ticks.x = element_line(), 
              axis.text = element_text(size = 11),
              # Facets
              strip.text = element_text(color = "white", 
                                        face = "bold",
                                        size = 11),
              strip.background = element_rect(
                color = "white", 
                fill = "#00736D", 
                size = 1, 
                linetype="solid"))
          )



# Building a GMRI theme based on Wall street Journal and NYTimes theme
# base settings from {ggthemes}
theme_gmri <- function(base_size = 10, 
                       bg_color = "lightblue", 
                       base_family = "sans", 
                       title_family = "sans",
                       facet_color = "teal") {
    # Color from gmRi palette, sets background color
    #colorhex <- gmRi::gmri_cols()[bg_color]
    facet_hex <- gmri_cols()[facet_color]
    
    # Set up theme
    theme_foundation(
      base_size = base_size, 
      base_family = base_family) + 
        theme(
          
          # Major Elements
          line = element_line(linetype = 1, colour = "black"), 
          rect = element_rect(fill = "transparent", 
                              linetype = 0, 
                              colour = NA), 
          text = element_text(colour = "black"), 
          title = element_text(family = title_family, size = 12), 
          
          # Axis elements
          axis.text.x = element_text(colour = NULL), 
          axis.text.y = element_text(colour = NULL), 
          axis.ticks = element_line(colour = NULL), 
          axis.ticks.y = element_blank(), 
          axis.ticks.x = element_line(colour = NULL), 
          axis.line = element_line(), 
          axis.line.y = element_blank(), 
          axis.text = element_text(size = 11),
          
          # Legend Elements
          legend.background = element_rect(), 
          legend.position = "top", 
          legend.direction = "horizontal", 
          legend.box = "vertical", 
          legend.title = element_text(size = 9),
          legend.text = element_text(size = 7.5),
          
          # Panel/Grid Setup
          panel.grid = element_line(colour = NULL, linetype = 3, color = "gray70"), 
          panel.grid.major = element_line(colour = "black"), 
          panel.grid.major.x = element_blank(), 
          panel.grid.minor = element_blank(), 
          
          # Title and Caption Details
          plot.title = element_text(hjust = 0, face = "bold", size = 14),
          plot.subtitle = element_text(size = 9),
          plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
          plot.margin = unit(c(1, 1, 1, 1), "lines"), 
          
          # Facet Details
          strip.text = element_text(color = "white", face = "bold", size = 11),
          strip.background = element_rect(
            color = "white", 
            fill = facet_hex, 
            size = 1, 
            linetype="solid"))
}

# Set the theme
# theme_set(theme_gmri(base_size = 10, base_family = "sans", title_family = "sans", facet_color = "teal"))

2021 Sea Surface Temperature

An Exceptional Year

2021 was an exceptionally warm year for the Gulf of Maine. It had the warmest fall on record for the area, and the second warmest summer. Much of the year the region experienced what would be considered “heatwave conditions”.

hide
# Load the timeseries
region_timeseries <- read_csv(timeseries_path, col_types = cols(), guess_max = 1e6)


# Clean up the data - add labels
region_timeseries <- region_timeseries %>% 
  distinct(time, .keep_all = T) %>% 
  mutate(time = as.Date(time),
         yr = year(time),
         month_num = month(time),
         month = month(time, label = T, abbr = T),
         week = lubridate::week(time),
         doy = yday(time),
         season = metR::season(month_num, lang = "en"),
         #Set up correct year for winters, they carry across into next year
         season_eng = case_when(
           season == "SON" ~ "Fall",
           season == "DJF" ~ "Winter",
           season == "MAM" ~ "Spring",
           season == "JJA" ~ "Summer"),
         season_yr = ifelse( (season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
         )

# Get heatwave statuses for each day:
# Uses area weighted sst by default
region_hw <- pull_heatwave_events(temperature_timeseries = region_timeseries)

# Format dates and seasons and heatwave outcomes
region_hw <- region_hw %>% 
  mutate(time = as.Date(time),
         yr = year(time),
         month_num = month(time),
         month = month(time, label = T, abbr = T),
         week = lubridate::week(time),
         doy = yday(time),
         season = metR::season(month_num, lang = "en"),
         season_eng = case_when(
           season == "SON" ~ "Fall",
           season == "DJF" ~ "Winter",
           season == "MAM" ~ "Spring",
           season == "JJA" ~ "Summer"),
         #Set up correct year for winters, they carry across into next year
         season_yr = ifelse(
           (season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
         )

# Set up heatwave status outcomes for fancy table features:
region_hw <- region_hw %>% 
  mutate(hw_outcomes = case_when(
    mhw_event == TRUE ~ 1,
    mcs_event == TRUE ~ 0,
    TRUE ~ 0.5))

Departure from Historic Conditions

hide
# Pathto ERSST on box
ersst_path <- cs_path("okn", "ersst")

# Load temps and anoms (1982-2011 climatology)
ersst_temp <- stack(str_c(ersst_path, "ERSSTv5_merged.nc"))
ersst_anom <- stack(str_c(ersst_path, "ERSSTv5_anom.nc"))

# Crop GoM
# function to mask
mask_nc <- function(ras_obj, mask_shape){
  
  # Get stats from ranks
  m1 <- raster::mask(rotate(ras_obj), mask_shape)
  m1 <- raster::crop(m1, mask_shape)
  return(m1)}

# Load the shape, mask ERSST with it
gom_shape <- read_sf(poly_path)
gom_ersst <- mask_nc(ras_obj = ersst_temp, mask_shape = gom_shape)
gom_ersst_anom <- mask_nc(ras_obj = ersst_anom, mask_shape = gom_shape)

#  Clean them up - temp
ersst_temp_tl <- as.data.frame(cellStats(gom_ersst, mean)) %>% 
  rownames_to_column(var = "Date") %>% 
  setNames(c("Date", "sst")) %>% 
  mutate(Date = str_sub(Date, 2, -1),
         Date = str_replace_all(Date, "[.]", "-"),
         Date = as.Date(Date), 
         yr = str_sub(Date, 1, 4))

# anomalies
ersst_anom_tl <- as.data.frame(cellStats(gom_ersst_anom, mean)) %>% 
  rownames_to_column(var = "Date") %>% 
  setNames(c("Date", "sst_anom")) %>% 
  mutate(Date = str_sub(Date, 2, -1),
         Date = str_replace_all(Date, "[.]", "-"),
         Date = as.Date(Date), 
         yr = str_sub(Date, 1, 4))


# Join them together
ersst_tl <- left_join(ersst_temp_tl, ersst_anom_tl) %>% 
  mutate(yr = as.numeric(yr))

# Get Yearly Average
ersst_yr <- ersst_tl %>% 
  group_by(yr) %>% 
  summarise(sst = mean(sst), 
            sst_anom = mean(sst_anom),
            .groups = "drop") %>% 
  mutate(yr = as.numeric(yr)) %>% 
  filter(yr <= 2019) 


# Same for OISST
oisst_yr <- region_hw %>% 
  filter(between(yr, 1982, 2021)) %>% 
  group_by(yr) %>% 
  summarise(sst = mean(sst), .groups = "drop")

# Make Splines
ersst_smooth <- as.data.frame(spline(ersst_yr$yr, ersst_yr$sst))
oisst_smooth <- as.data.frame(spline(oisst_yr$yr, oisst_yr$sst))


# Plot them
ggplot() +
  geom_line(data = ersst_smooth, aes(x, y, color = "ERSSTv5"),
            group = 1, alpha = 0.4, linetype = 1) +
  geom_point(data = ersst_yr, aes(yr, sst, color = "ERSSTv5"), 
             size = 1) +
 geom_line(data = oisst_smooth, aes(x, y, color = "OISSTv2"),
            group = 1, alpha = 0.4, linetype = 1) +
  geom_point(data = oisst_yr, aes(yr, sst, color = "OISSTv2"), 
             size = 1) +
  scale_color_gmri() +
  theme_gmri() +
  labs(x = "", y = "Sea Surface Temperature", 
       color = "Temperature Data Record",
       title = "Long-Term Temperature Record for the Gulf of Maine",
       subtitle = "Current temperatures above 150-year highs") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
  facet_zoom(xlim = c(1982, 2021)) +
  theme(legend.position = "bottom")

Long term temperature record for the Gulf of Maine

Relationship to Climate Patterns

To better get a sense of what kind of year 2021 was it is helpful to relate it to what we would expect from similar years. One large factor in global weather patterns is the El Niño Southern Oscillation (ENSO).

hide
# Source: https://www.ncdc.noaa.gov/teleconnections/enso/soi


# library(rsoi)
# clim_idx <- download_soi()
soi_std <- tibble::tribble(#####
~"YEAR", ~"JAN", ~"FEB", ~"MAR", ~"APR", ~"MAY", ~"JUN", ~"JUL", ~"AUG", ~"SEP", ~"OCT", ~"NOV", ~"DEC",    
1951,   1.5,   0.9,  -0.1,  -0.3,  -0.7,   0.2,  -1.0,  -0.2,  -1.1,  -1.0,  -0.8,  -0.7,
1952,  -0.9,  -0.6,   0.5,  -0.2,   0.8,   0.7,   0.5,   0.1,  -0.2,   0.4,   0.0,  -1.2,
1953,   0.3,  -0.5,  -0.2,   0.2,  -1.7,   0.1,  -0.0,  -1.2,  -1.2,   0.1,  -0.3,  -0.5,
1954,   0.7,  -0.3,   0.3,   0.6,   0.5,   0.1,   0.4,   1.1,   0.2,   0.3,   0.1,   1.4,
1955,  -0.5,   1.9,   0.6,  -0.1,   1.0,   1.3,   1.6,   1.5,   1.3,   1.5,   1.2,   1.0,
1956,   1.3,   1.6,   1.3,   0.9,   1.4,   1.1,   1.1,   1.2,   0.1,   1.8,   0.2,   1.1,
1957,   0.6,  -0.1,   0.2,   0.2,  -0.7,   0.2,   0.2,  -0.5,  -0.9,   0.0,  -1.0,  -0.3,
1958,  -1.9,  -0.5,   0.3,   0.4,  -0.5,   0.3,   0.4,   0.9,  -0.3,   0.1,  -0.4,  -0.6,
1959,  -0.9,  -1.4,   1.3,   0.4,   0.5,  -0.1,  -0.3,  -0.1,   0.0,   0.5,   0.9,   0.9,
1960,   0.1,   0.1,   1.0,   0.8,   0.5,   0.1,   0.5,   0.8,   0.7,   0.1,   0.5,   0.8,
1961,  -0.3,   0.9,  -1.8,   0.8,   0.3,   0.1,   0.2,   0.2,   0.1,  -0.3,   0.5,   1.5,
1962,   2.0,  -0.3,   0.1,   0.2,   1.1,   0.7,   0.1,   0.6,   0.4,   1.0,   0.3,   0.2,
1963,   1.0,   0.6,   1.1,   0.8,   0.4,  -0.5,  -0.1,   0.0,  -0.6,  -1.2,  -0.8,  -1.2,
1964,  -0.4,   0.0,   1.1,   1.1,   0.2,   0.8,   0.6,   1.5,   1.3,   1.3,   0.2,  -0.3,
1965,  -0.4,   0.4,   0.8,  -0.5,   0.2,  -0.6,  -1.8,  -0.7,  -1.3,  -0.9,  -1.5,   0.2,
1966,  -1.3,  -0.2,  -0.9,  -0.2,  -0.4,   0.3,   0.1,   0.6,  -0.2,  -0.1,  -0.0,  -0.3,
1967,   1.7,   1.7,   1.2,  -0.0,  -0.0,   0.6,   0.2,   0.7,   0.5,   0.1,  -0.4,  -0.6,
1968,   0.5,   1.3,   0.1,   0.0,   1.2,   1.1,   0.7,   0.3,  -0.3,  -0.1,  -0.3,   0.2,
1969,  -1.5,  -0.5,   0.4,  -0.4,  -0.2,   0.2,  -0.5,  -0.1,  -1.0,  -0.9,  -0.1,   0.4,
1970,  -1.1,  -1.0,   0.6,  -0.1,   0.4,   1.0,  -0.4,   0.6,   1.2,   1.0,   1.6,   1.9,
1971,   0.4,   2.0,   2.3,   1.7,   0.9,   0.4,   0.2,   1.5,   1.4,   1.7,   0.5,   0.3,
1972,   0.5,   1.1,   0.6,  -0.1,  -1.6,  -0.5,  -1.4,  -0.5,  -1.4,  -0.9,  -0.3,  -1.3,
1973,  -0.3,  -1.4,   0.7,   0.1,   0.4,   1.1,   0.6,   1.3,   1.2,   0.8,   2.6,   1.8,
1974,   2.4,   2.1,   2.4,   0.9,   1.0,   0.4,   1.1,   0.8,   1.1,   0.9,  -0.1,   0.2,
1975,  -0.5,   0.8,   1.6,   1.2,   0.6,   1.3,   1.9,   2.0,   2.1,   1.7,   1.2,   2.1,
1976,   1.4,   1.7,   1.7,   0.3,   0.4,   0.3,  -0.9,  -0.8,  -1.1,   0.4,   0.7,  -0.3,
1977,  -0.4,   1.2,  -0.5,  -0.4,  -0.5,  -0.9,  -1.1,  -0.8,  -0.8,  -1.0,  -1.3,  -1.1,
1978,  -0.3,  -2.7,  -0.2,  -0.3,   1.4,   0.7,   0.6,   0.4,   0.1,  -0.4,  -0.0,  -0.1,
1979,  -0.4,   1.0,   0.1,  -0.1,   0.5,   0.6,   1.3,  -0.2,   0.1,  -0.1,  -0.4,  -0.7,
1980,   0.4,   0.3,  -0.4,  -0.6,  -0.0,  -0.0,  -0.0,   0.4,  -0.5,   0.0,  -0.3,  -0.1,
1981,   0.4,  -0.2,  -1.3,  -0.1,   0.8,   1.2,   0.8,   0.7,   0.3,  -0.4,   0.2,   0.5,
1982,   1.2,   0.3,   0.6,   0.1,  -0.3,  -1.0,  -1.5,  -1.7,  -1.7,  -1.7,  -2.6,  -2.2,
1983,  -3.5,  -3.6,  -2.4,  -0.9,   0.6,   0.0,  -0.6,   0.1,   0.9,   0.4,  -0.1,   0.0,
1984,   0.2,   0.9,  -0.2,   0.3,   0.2,  -0.3,   0.2,   0.4,   0.1,  -0.3,   0.3,  -0.1,
1985,  -0.3,   1.2,   0.8,   1.2,   0.4,  -0.4,  -0.1,   1.0,   0.0,  -0.4,  -0.2,   0.2,
1986,   1.0,  -1.0,   0.5,   0.3,  -0.2,   1.0,   0.3,  -0.4,  -0.5,   0.6,  -1.2,  -1.4,
1987,  -0.7,  -1.2,  -1.3,  -1.4,  -1.3,  -1.1,  -1.4,  -0.9,  -1.0,  -0.4,  -0.0,  -0.5,
1988,  -0.1,  -0.4,   0.6,   0.1,   0.9,   0.1,   1.0,   1.5,   1.8,   1.4,   1.7,   1.2,
1989,   1.5,   1.2,   1.1,   1.6,   1.2,   0.7,   0.9,  -0.3,   0.5,   0.8,  -0.2,  -0.5,
1990,  -0.1,  -1.8,  -0.4,   0.2,   1.2,   0.3,   0.5,  -0.2,  -0.7,   0.3,  -0.5,  -0.2,
1991,   0.6,   0.3,  -0.7,  -0.6,  -1.0,  -0.1,   0.0,  -0.4,  -1.5,  -1.0,  -0.7,  -1.8,
1992,  -2.9,  -0.9,  -2.0,  -1.0,   0.3,  -0.6,  -0.6,   0.4,   0.1,  -1.4,  -0.7,  -0.6,
1993,  -0.9,  -0.7,  -0.5,  -1.2,  -0.3,  -0.8,  -0.8,  -0.9,  -0.7,  -1.1,  -0.1,   0.2,
1994,  -0.1,   0.3,  -0.7,  -1.3,  -0.7,  -0.4,  -1.3,  -1.2,  -1.6,  -1.1,  -0.6,  -1.2,
1995,  -0.4,  -0.1,   0.8,  -0.7,  -0.4,   0.1,   0.4,   0.3,   0.3,   0.0,   0.0,  -0.5,
1996,   1.0,   0.3,   1.1,   0.8,   0.3,   1.2,   0.7,   0.7,   0.6,   0.6,  -0.1,   0.9,
1997,   0.5,   1.7,  -0.4,  -0.6,  -1.3,  -1.4,  -0.8,  -1.4,  -1.4,  -1.5,  -1.2,  -1.0,
1998,  -2.7,  -2.0,  -2.4,  -1.4,   0.3,   1.0,   1.2,   1.2,   1.0,   1.1,   1.0,   1.4,
1999,   1.8,   1.0,   1.3,   1.4,   0.2,   0.3,   0.5,   0.4,  -0.1,   1.0,   1.0,   1.4,
2000,   0.7,   1.7,   1.3,   1.2,   0.4,  -0.2,  -0.2,   0.7,   0.9,   1.1,   1.8,   0.8,
2001,   1.0,   1.7,   0.9,   0.2,  -0.5,   0.3,  -0.2,  -0.4,   0.2,  -0.0,   0.7,  -0.8,
2002,   0.4,   1.1,  -0.2,  -0.1,  -0.8,  -0.2,  -0.5,  -1.0,  -0.6,  -0.4,  -0.5,  -1.1,
2003,  -0.2,  -0.7,  -0.3,  -0.1,  -0.3,  -0.6,   0.3,   0.1,  -0.1,   0.0,  -0.3,   1.1,
2004,  -1.3,   1.2,   0.4,  -0.9,   1.0,  -0.8,  -0.5,  -0.3,  -0.3,  -0.1,  -0.7,  -0.8,
2005,   0.3,  -3.1,   0.3,  -0.6,  -0.8,   0.4,   0.2,  -0.3,   0.4,   1.2,  -0.2,  -0.0,
2006,   1.7,   0.1,   1.8,   1.1,  -0.5,  -0.2,  -0.6,  -1.0,  -0.6,  -1.3,   0.1,  -0.3,
2007,  -0.8,  -0.1,   0.2,  -0.1,  -0.1,   0.5,  -0.3,   0.4,   0.2,   0.7,   0.9,   1.7,
2008,   1.8,   2.6,   1.4,   0.7,  -0.1,   0.6,   0.3,   1.0,   1.2,   1.3,   1.3,   1.4,
2009,   1.1,   1.9,   0.4,   0.8,  -0.1,   0.1,   0.2,  -0.2,   0.3,  -1.2,  -0.6,  -0.7,
2010,  -1.1,  -1.5,  -0.7,   1.2,   0.9,   0.4,   1.8,   1.8,   2.2,   1.7,   1.3,   2.9,
2011,   2.3,   2.7,   2.5,   1.9,   0.4,   0.2,   1.0,   0.4,   1.0,   0.8,   1.1,   2.5,
2012,   1.1,   0.5,   0.7,  -0.3,   0.0,  -0.4,  -0.0,  -0.2,   0.2,   0.3,   0.3,  -0.6,
2013,  -0.1,  -0.2,   1.5,   0.2,   0.8,   1.2,   0.8,   0.2,   0.3,  -0.1,   0.7,   0.1,
2014,   1.4,   0.1,  -0.9,   0.8,   0.5,   0.2,  -0.2,  -0.7,  -0.7,  -0.6,  -0.9,  -0.6,
2015,  -0.8,   0.2,  -0.7,  -0.0,  -0.7,  -0.6,  -1.1,  -1.4,  -1.6,  -1.7,  -0.5,  -0.6,
2016,  -2.2,  -2.0,  -0.1,  -1.2,   0.4,   0.6,   0.4,   0.7,   1.2,  -0.3,  -0.1,   0.3,
2017,   0.2,  -0.1,   0.9,  -0.2,   0.3,  -0.4,   0.8,   0.5,   0.6,   0.9,   0.9,  -0.1,
2018,   1.1,  -0.5,   1.5,   0.5,   0.4,  -0.1,   0.2,  -0.3,  -0.9,   0.4,  -0.1,   1.0,
2019,  -0.0,  -1.4,  -0.3,   0.1,  -0.4,  -0.5,  -0.4,  -0.1,  -1.2,  -0.4,  -0.8,  -0.6,
2020,   0.2,  -0.1,  -0.1,   0.2,   0.4,  -0.4,   0.4,   1.1,   0.9,   0.5,   0.7,   1.8,
2021,   1.9,   1.5,   0.4,   0.3,   0.5,   0.4,   1.4,   0.6,   0.8,   0.7,   1.0,   1.5
  
)

#####

# Pick a date centered in each month so you can plot them as one timeline
month_centers <- data.frame(Month = toupper(month.abb),
                            Month_num = str_pad(seq(1,12,1), width = 2, pad = "0", side = "left"))
                            
# Reformat
soi_long <- soi_std %>% 
  pivot_longer(names_to = "Month", values_to = "soi_z", cols = JAN:DEC) %>% 
  left_join(month_centers, by = "Month") %>% 
  mutate(Date = as.Date(str_c(YEAR, "-", Month_num, "-15")),
         ribbon_fill = ifelse(soi_z > 0,"above", "below"),
         ribbon_above = ifelse(soi_z > 0, soi_z, NA),
         above_ymin   = ifelse(soi_z > 0, 0, NA),
         ribbon_below = ifelse(soi_z < 0, soi_z, NA),
         below_ymin = ifelse(soi_z < 0, 0, NA)
         )

# The plotalllayers thing wasnt cutting it
# Plot the ribbon plot instead
cutpoints <- rev(seq(-4, 4, 1))
col_rdbu <- rev(RColorBrewer::brewer.pal(n = length(cutpoints), name = "RdBu"))


# fix GeomRibbon
GeomRibbon$handle_na <- function(data, params) {  data }

# Make figure
ggplot(filter(soi_long, YEAR >= 1982), aes(x = Date, y =  soi_z)) +
  geom_ribbon(aes(ymin = above_ymin, ymax = ribbon_above, fill = "La Niña"),  alpha = 0.8) +
  geom_ribbon(aes(ymin = below_ymin, ymax = ribbon_below, fill = "El Niño"), alpha = 0.8) +
  scale_fill_manual(values = c("La Niña" = col_rdbu[2], "El Niño" = col_rdbu[7])) +
  geom_line(size = 0.25, color = "gray20") + 
  labs(x = "Date", y = "SOI", title = "ENSO Index",
       caption = "Years overlapping OISST data displayed", fill = "") + 
  theme_gmri() +
  theme(legend.position = "bottom")

When we look at average temperatures for the Gulf of Maine, and how that relates to the ENSO index, we found there was NO relationship. Meaning that a current month’s El niño/ La Nińa status did not correlate with warmer/hotter months.

A similar climate index to ENSO, with more relevance to our area is the North Atlantic Oscillation (NAO).

Looking at the NAO index, we found there was ___ relationship. Meaning that a current month’s El niño/ La Nińa status did ____ with warmer/hotter months.

Daily Temperature Record

hide
 # Last Calendar Year
this_yr <- region_hw %>% 
  filter(year(time) == 2021)


# Set colors by name
color_vals <- c(
  "Sea Surface Temperature" = "royalblue",
  "Heatwave Event"          = "darkred",
  "MHW Threshold"           = "coral3",
  "Daily Climatology"        = "gray30")


# Set the label with degree symbol
ylab <- "Sea Surface Temperature"



# Plot the last 365 days
hw_temp_p <- ggplot(this_yr, aes(x = time)) +
    geom_segment(aes(x = time, xend = time, y = seas, yend = sst), 
                 color = "royalblue", alpha = 0.25) +
    geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe), 
                 color = "darkred", alpha = 0.25) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    # geom_textpath(aes(y = mhw_thresh,  color = "MHW Threshold"), label = "Heatwave Threshold", hjust = 0.8, lty = 3) +
    #geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
    geom_textpath(aes(y = seas,  color = "Daily Climatology"), label = "Climatology Mean", hjust = 0.5, lty = 2) +
    scale_color_manual(values = color_vals) +
    scale_x_date(date_labels = "%b", date_breaks = "1 month") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
                                         labels =  number_format(suffix = " \u00b0F")),
                     labels = number_format(suffix = " \u00b0C")) +
    theme(legend.title = element_blank(),
          legend.position = "none") +
    labs(x = "", 
         y = ylab)



# Plot the last 365 days - anomaly scale
linetype_key <- c(
  "Sea Surface Temperature Anomaly" = 1,
  "Heatwave Event"                  = 1,
  "MHW Threshold"                   = 3,
  "Daily Climatology"               = 2)

# Same Plot as Anomalies
hw_anom_p <- this_yr %>% 
  mutate(
    sst = sst,
    seas = seas,
    sst_anom = sst_anom,
    mhw_thresh = mhw_thresh,
    anom_thresh = mhw_thresh - seas,
    anom_hwe = hwe - seas) %>% 
ggplot(aes(x = time)) +
    geom_segment(aes(x = time, xend = time, 
                     y = 0, yend = sst_anom), 
                 color = "royalblue", alpha = 0.25) +
    geom_segment(aes(x = time, xend = time, 
                       y = anom_thresh, yend = anom_hwe), 
                 color = "darkred", alpha = 0.25) +
    geom_line(aes(y = sst_anom, color = "Sea Surface Temperature Anomaly")) +
    geom_line(aes(y = anom_hwe, color = "Heatwave Event")) +
    geom_line(aes(y = anom_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    #geom_line(aes(y = 0, color = "Daily Climatology"), lty = 2, size = .5) +
    geom_textpath(aes(y = 0,  color = "Daily Climatology"), 
                  label = "Climatology Mean", hjust = 0.5, lty = 2, show.legend = FALSE) +
    scale_color_manual(values = color_vals) +
    scale_linetype_manual(values = linetype_key, guide = "none") +
    scale_x_date(date_labels = "%b", date_breaks = "1 month") +
    guides(color = guide_legend(override.aes = list(linetype = linetype_key), nrow = 2)) +
    theme(legend.title = element_blank(),
          legend.position = "top") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
                                           labels =  number_format(suffix = " \u00b0F")),
                       labels = number_format(suffix = " \u00b0C")) +
    labs(x = "", 
         y = "Temperature Anomaly", 
         caption = paste0("Climate reference period :  1982-2011"))


# Show figure
hw_temp_p / hw_anom_p

Timelines of observed temperature and observed temperature anomalies for the year of 2021, and how they compare to the mean and 90th percentile of the climate reference period

2021 Heatwave Events

hide
# Isolate 2021
hw_2021 <- filter(region_hw, yr == "2021")

# Get Event numbers
events_21 <-  hw_2021 %>% drop_na(mhw_event_no) %>% distinct(mhw_event_no) %>% pull(mhw_event_no)
events_num <-  length(events_21)


# Pull their data
events_dat <- filter(region_hw, mhw_event_no %in% as.character(events_21))
  

# How long on average
avg_duration <- round(nrow(events_dat) / events_num, digits = 0)

During 2021 there were three main marine heatwave events, including one that began on the new year. The average duration for each heatwave event was 75 days, or over three months long.

Heatwave Seasonality

hide
# Set new axis dimensions, y = year, x = day within year
# use a flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_hw %>% 
  mutate(year = year(time),
         yday = yday(time),
         flat_date = as.Date(yday-1, origin = base_date))


# Polar plot comparing when each year experienced peak heatwave conditions
polar_dat <- grid_data %>% filter(year %in% c("2012", "2021")) %>% 
  mutate(hw_line =  ifelse(mhw_event == TRUE, sst_anom, NA),
         hw_mag = ifelse(year == "2012", hw_line * -1, hw_line))

# Ribbon?
polar_dat  %>%
ggplot(aes(flat_date, y = hw_line, ymin = 0, ymax = hw_line)) +
  geom_ribbon(fill = "darkred", alpha = 0.75) +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) + 
  facet_wrap( ~ year) + 
  geomtextpath::coord_curvedpolar() +
  # theme_gmri() +
  labs(x = "", y = "Temperature Anomalies")

hide
# Side by side - need data separately
polar_dat %>% 
  ggplot(aes(x = hw_line, xend = 0, 
             y = flat_date, yend = flat_date, color = factor(year))) +
  geom_segment(alpha = 0.4) +
  scale_y_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0), ) + 
  facet_wrap( ~ year) + 
  scale_color_gmri() +
  theme_gmri() +
  # coord_flip() +
  labs(x = "", y = "Temperature Anomalies", fill = "Year")

Following Each Heatwave

Do some breakdowns of how each heatwave looked: - where was heat concentrated - how long did it last - was anything unique about it?

Heatwave Frequency Heatmap

hide
# Prep the legend title
guide_lab <- "Sea Surface Temperature Anomaly \u00b0C"

# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(grid_data$sst_anom)) * c(-1, 1)


# Assemble heatmap plot
heatwave_heatmap <- ggplot(grid_data, aes(x = flat_date, y = year)) +
  
  # background box fill for missing dates
  geom_rect(xmin = base_date, xmax = base_date + 365, 
            ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5, 
            fill = "gray75", color = "transparent") +
  
  # tile for sst colors
  geom_tile(aes(fill = sst_anom)) +
  
  # points for heatwave events
  geom_point(data = filter(grid_data, mhw_event == TRUE),
             aes(x = flat_date, y = year), size = .25)  +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
  scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
  labs(x = "", 
       y = "",
       "\nClimate reference period : 1982-2011") +
  
  #scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", 
                       limit = limit) +
  
  #5 inches is default rmarkdown height for barheight
  guides("fill" = guide_colorbar(title = guide_lab, 
                                 title.position = "right", 
                                 title.hjust = 0.5,
                                 barheight = unit(4.8, "inches"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +  
  theme_classic() +
  theme(legend.title = element_text(angle = 90))


# Assemble pieces
heatwave_heatmap

Heatmap of heatwave temperatures and heatwave days

Temperature Horizons

hide
# make slight modification to grid_data
hori_data <- grid_data %>% 
  mutate(year = factor(year),
         year = fct_rev(year))


# cut out outliers and set cutpoints
cutpoints <- hori_data %>% 
  mutate(
    outlier = between(
      sst_anom,
      quantile(sst_anom, 0.25, na.rm=T)-
        1.5*IQR(sst_anom, na.rm=T),
      quantile(sst_anom, 0.75, na.rm=T)+
        1.5*IQR(sst_anom, na.rm=T))) %>% 
  filter(outlier)

# Set origin
ori <- 0

# Manually pick cut points
sca <- seq(-4, 4, length.out = 9)

sca_bins <- str_c(sca[1:(length(sca)-1)], " to ", sca[2:length(sca)], "\u00b0C")
sca_labels <- rev(sca_bins)


# Function to plot one or more years
anom_horizon_plot <- function(hori_data, origin = ori, scale_cutpoints = sca, labels = sca_labels){
  
   ggplot(hori_data) +
  geom_horizon(aes(flat_date, 
                   sst_anom,
                   fill = ..Cutpoints..), 
               origin = ori, 
               horizonscale = sca) +
  scale_fill_hcl(palette = 'RdBu', reverse = T, labels = sca_labels) +
  facet_grid(year~.) +
  theme(
    panel.grid = element_blank(),
    panel.spacing.y=unit(0, "lines"),
    strip.text.y = element_text(size = 7, angle = 0, hjust = 0),
    axis.text.y = element_blank(),
    axis.title.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.border = element_blank(),
    legend.position = "left"
    ) +
  scale_x_date(expand=c(0,0), 
               date_breaks = "1 month", 
               date_labels = "%b") +
  labs(x = '', 
       fill = "Temperature Anomaly \u00b0C") 
  
  
}


# Plot Horizon Plot Using All Years
anom_horizon_plot(hori_data = filter(hori_data, yr <= 2021), 
                  origin = ori, 
                  scale_cutpoints = sca) +
  labs(title = "Temperature Anomaly Horizons",
       subtitle = "How deviation from the norm has become the new-normal",
       caption = "Climatatological Reference Period: 1982-2011")

Horizon plot of all years and their temperature anomaly horizons

Comparing the Gulf’s Two Hottest Years

Early into 2021 it was apparent that the year was on-par with the previous title-holder for warmest year on record.

hide
# Reset the outliers
# cut out outliers and set cutpoints
cutpoints_2 <- hori_data %>% 
  filter(year %in% c("2012", "2021")) %>% 
  mutate(
    outlier = between(
      sst_anom,
      quantile(sst_anom, 0.25, na.rm=T)-
        1.5*IQR(sst_anom, na.rm=T),
      quantile(sst_anom, 0.75, na.rm=T)+
        1.5*IQR(sst_anom, na.rm=T))) %>% 
  filter(outlier)



# 2012
hori_2012 <- hori_data %>% 
  filter(year == "2012") %>% 
  anom_horizon_plot(scale_cutpoints = cutpoints_2) +
  facet_wrap(~year, strip.position = "top")
  
# 2021
hori_2021 <- hori_data %>% filter(year == "2021") %>% 
  anom_horizon_plot(scale_cutpoints = cutpoints_2) +
  facet_wrap(~year, strip.position = "top") +
  theme(plot.caption = element_text(colour = "gray25", size = 6)) +
  labs(caption = "Climatatological Reference Period: 1982-2011")

# Put them together as one figure
comparison_horizons <- hori_2012 / hori_2021 + plot_layout(guides = "collect")

# Do a bunch of formatting
comparison_horizons &
  theme(
    #Legend
    legend.position = "top",
    legend.key.height = unit(.5, "lines"),
    legend.key.width = unit(3.75, "lines")) &
   guides(fill = guide_legend(title = "Temperature Anomaly",
                              title.hjust = 0.5, 
                              title.position = "top",
                              label.position = "bottom", 
                              nrow = 1, 
                              byrow = T, reverse = T)) &
  plot_annotation(title = "Comparing the Gulf of Maine's Two Hottest Years", 
                  subtitle = "Duration of temperature anomaly horizons, colored by their strength")
  grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

Single year horizon plots comparing 2012 and 2022

Streamplot

hide
# Set Factor Order for Bins
bin_order <- c(
  "Less than -1", 
  "-1 to -0.5",
  "Within 0.5",
  "0.5 to 1", 
  "Greater than 1")

# Set Colors
bin_colors <- RColorBrewer::brewer.pal(n = length(bin_order), 
                                       name = "RdBu")
bin_colors <- rev(bin_colors)

# Count Days in Bins
day_types <- region_hw %>% 
  mutate(
    `Less than -1`   = if_else(sst_anom <= -1, TRUE, FALSE),
    `-1 to -0.5`     = if_else(between(sst_anom, -1, -0.5), TRUE, FALSE),
    `Within 0.5`     = if_else(between(sst_anom, -0.5, 0.5), TRUE, FALSE),
    `0.5 to 1`       = if_else(between(sst_anom, 0.5, 1), TRUE, FALSE),
    `Greater than 1` = if_else(sst_anom >= 1, TRUE, FALSE)) %>% 
  group_by(yr) %>% 
  summarise(
    `Less than -1`   = sum(`Less than -1`, na.rm = T),
     `-1 to -0.5`    = sum(`-1 to -0.5`, na.rm = T),
    `Within 0.5`     = sum(`Within 0.5`, na.rm = T),
    `0.5 to 1`       = sum(`0.5 to 1`, na.rm = T),
    `Greater than 1` = sum(`Greater than 1`, na.rm = T),
    .groups = "drop") %>% 
  pivot_longer(cols = "Less than -1":"Greater than 1", names_to = "Anomaly Strength", values_to = "day_total", values_drop_na = FALSE) %>% 
  mutate(day_total = if_else(is.na(day_total), 0L, day_total),
         `Anomaly Strength` = factor(`Anomaly Strength`, levels = bin_order))


# Plot the Streamplot
hw_streamplot <- ggplot(day_types, aes(yr, y = day_total, fill = `Anomaly Strength`)) +
  geom_stream(type = "proportion") +
  #geom_segment(y = .5, yend = 0.5, linetype = 3, color = "gray50", x = 1981, xend = 2021, size = 0.25) +
  scale_x_continuous(breaks = seq(1980,2030, by = 10), 
                     minor_breaks = seq(1975, 2030, by = 5), 
                     expand = c(.03,0)) +
  scale_y_continuous(labels = scales::percent_format(), expand = c(0,0)) + 
  scale_fill_manual(values = setNames(bin_colors, bin_order)) + 
  labs(y = "", x = "", fill = "", 
       title = "A Warmer Gulf of Maine",
       subtitle = "Fractions of Each Year Experienced as Hot/Cold Anomalies",
       caption = "1981 not a complete years in the data*") +
  theme_minimal(base_size = 10) +
  theme(
    axis.line.x = element_line(),
    axis.ticks.x = element_line(),
    axis.line.y = element_line(),
    axis.ticks.y = element_line(),
    legend.position = "top",
    legend.key.height = unit(.5, "lines"),
    legend.key.width = unit(3.75, "lines"),
    panel.grid.minor = element_blank(),
    panel.grid = element_line(size = 0.2),
    plot.margin = unit(c(.5,1,.3,.5), "cm"),
    plot.title.position = "plot",
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 9),
    plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
    legend.title = element_text(size = 9),
    legend.text = element_text(size = 7.5),
    legend.margin = margin(c(7,0,0,0)),
    legend.justification = "center") +
   guides(fill = guide_legend(title = "Anomaly Strength \u00b0C",
                             title.hjust = 0.5, 
                             title.position = "left",
                             label.position = "bottom", 
                             nrow = 1, 
                             byrow = T))


# Show figure
hw_streamplot
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))

Streamplot tracking fraction of each year spent in varying degrees of temperature anomalies

Characteristics of the Hottest Years

Multivariate figure comparing characteristics of both years: 2012 & 2021

hide
# Use all years so they can be ranked independently

# Metrics?
year_mets <- hori_data %>% 
  filter(year != "2022") %>% 
  group_by(year) %>% 
  summarise(
    total_days                         = n(),
    `Average Temperature`              = round(mean(sst, na.rm = T), 2),
    `Temperature Above Average`        = round(mean(sst_anom, na.rm = T), 2),
    `Cumulative Degrees Above Average` = round(sum(sst_anom, na.rm = T), 0),
    hw_days                            = sum(mhw_event, na.rm = T),
    hw_events                          = n_distinct(mhw_event_no),
    peak_temp                          = max(sst),
    peak_anom                          = max(sst_anom)) %>% 
  mutate(
    hw_day_pct = round((hw_days/total_days) * 100, 2),
    `Average Heatwave Length` = round(hw_days/hw_events,  1)
  ) %>% 
  pivot_longer(names_to = "Var", values_to = "Metric", cols = `Average Temperature`:`Average Heatwave Length`)


# Snipe out the top5 for each
year_met_ranks <- year_mets %>% 
  filter(Var %in% c("Average Temperature", "Temperature Above Average", "Average Heatwave Length", "Cumulative Degrees Above Average")) %>% 
  group_by(Var) %>% 
  slice_max(n = 5, order_by = Metric) %>% 
  ungroup() %>%
  mutate(
    yr_focus = ifelse(year == "2021", TRUE, FALSE),
    year = reorder_within(year, Metric, Var)) 


year_met_ranks %>% 
  ggplot(aes(year, Metric, color = yr_focus)) +
    geom_col(aes(fill = yr_focus), show.legend = F) +
    geom_label(aes(label = Metric),  show.legend = FALSE) +
    facet_wrap(~Var, scales = "free", ncol = 2) +
    coord_flip() +
    scale_x_reordered() +
    scale_y_continuous(expand = expansion(mult = c(0, .2))) +
    labs(y = "",
         x = NULL,
         title = "How 2021 Stacks Up:",
         subtitle = "Ranking Characteristics of Acute and Prolonged Warming Events") +
  theme(plot.title = element_text(hjust = 0)) + 
  scale_color_gmri() +
  scale_fill_gmri() +
  theme(panel.border = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major. = element_line(linetype = 3, color = "gray50"))

4-panel plot comparing different features of hot years ranked by different metrics

Maps

Biological Events

 

A work by Adam A. Kemberling

Akemberling@gmri.org